;+
; NAME:
;   PROB_GRSMATCH
; PURPOSE:
;   Generate a spectrum for sources in the GRS overlap by integrating over
;   source and subracting off a boundary spectrum.
; CATEGORY:
;   distance-omnibus
; CALLING SEQUENCE:
;   probability = PROB_GRSMATCH(bdist, [bgps = bgps, grs_dir = grs_dir, $
;      bgps_dir = bgps_dir, map_dir = map_dir, obj_dir = obj_dir, spectrum =
;      spectrum)  
;
;
; INPUTS:
;   BDIST -- BGPS distance structure
;            BGPS -- Original BGPS structure to use.  Defaults to
;                    'bgps.v11.sav' in the BGPS_DIR/catalogs directory
; KEYWORD PARAMETERS:
;   GRS_DIR -- Directory containing GRS spectral line cubes.
;   BGPS_DIR -- Directory containing the BGPS data
;   MAP_DIR -- Directory containing the BGPS IMAGE files
;   OBJ_DIR -- Directory containing the BGPS object files
; OUTPUTS:
;   PROBABILITY -- probability vector normalized to unity.
; OPTIONAL OUTPUTS:
;   SPECTRUM -- Integrated masked spectrum used for probability calculation.
; MODIFICATION HISTORY:
;
;       Thu Jul 22 11:14:30 2010, erosolo <erosolo@>
;
;		Updated to newest version.
;
;       Fri Jun 11 13:36:03 2010, erosolo <erosolo@>
;                 Documented
;-

function prob_grsmatch, dstr, bgps = bgps, spectrum = spectrum,$
                        grs_dir = grs_dir, bgps_dir = bgps_dir,$
                        obj_dir = obj_dir, map_dir = map_dir

  COMPILE_OPT OBSOLETE
  
  if n_elements(grs_dir) eq 0 then grs_dir = './local/grs/'
  if n_elements(bgps_dir) eq 0 then bgps_dir = './local/bgps/'
  if n_elements(obj_dir) eq 0 then obj_dir = bgps_dir+'label/'
  if n_elements(map_dir) eq 0 then map_dir = bgps_dir+'maps/'


  defsysv, '!MW', exists = exists
  if not exists then galactic_params 

  R0 = !mw.r0
  V0 = !mw.v0

  vs = 15.0    ;; Correction for HMSFR
  if n_elements(errvlsr) eq 0 then errvlsr = 7.0 
; Assume a 7 km/s error in the velocities

  if n_elements(dvec) ne 0 then  d = dvec else $
     d = dindgen(!mw.nbins)*!mw.binsize+!mw.binstart
  prob = dblarr(!mw.nbins)+1d0/!mw.nbins

  if n_elements(bgps) eq 0 then $ 
     restore,file=bgps_dir+'catalogs/bgps.v11.sav'


  hit = where(dstr.cnum eq bgps.cloudnum,ct)
  if ct ne 1 then return, prob
  bgps = bgps[hit]


; Forward compatibility for ammonia spectral structures
  flag = total(tag_names(bgps) eq 'BGPS_FILENAME')
  if flag then filename = bgps.bgps_filename else $
     filename = bgps.filename
  catroots = (filename)
  slashpos = strpos(catroots, '/', /reverse_search)
  catroots = strmid(catroots, slashpos[0]+1, 60)
  
; Initialize structuring element
  r = 2
  elt = shift(dist(2*r+1,2*r+1),r,r) le r

  lastfn = ''
  lastbgps = ''

  if bgps.glon lt 14 or bgps.glon gt 60 then return, prob

  if lastbgps ne catroots then begin
     data = readfits(map_dir+catroots, bolo_hd)
     extast, bolo_hd, bolo_astrom
     lastbgps = catroots
     sz = size(data)
     char = stregex(catroots, '_map')
     labelname = strmid(catroots, 0, char+1)+$
                 'label'+strmid(catroots, char+1, 30)
     obj = readfits(obj_dir+labelname+'.gz')
  endif

; Determine where BOLOCAM points land in GRS data

;     ad2xy,bgps.glon_max,bgps.glat_max,bolo_astrom,xpix,ypix;
;     label = obj[xpix,ypix]
  label = bgps.cloudnum
  objmask = obj eq label

  ind = where(objmask, ct)
  if ct eq 0 then return, prob
  wts = data[ind]
  x_orig = ind mod sz[1]
  y_orig = ind / sz[1]

  border = (dilate(objmask, elt)-objmask)*(obj eq 0)
  ind2 = where(border,ctbdr)
  x_bdr = ind2 mod sz[1]
  y_bdr = ind2 / sz[1]
  
  xy2ad, x_orig, y_orig, bolo_astrom, glon, glat
  xy2ad, x_bdr, y_bdr, bolo_astrom, glonbdr, glatbdr

  wtorder = sort(glon)
  glon = glon[wtorder]
  glat = glat[wtorder]
  wts = wts[wtorder]
  lout = strcompress(string(round(glon)), /rem)
  filename = 'grs-'+lout+'-cube.fits'
  fn = grs_dir+filename

  wt2order = sort(glonbdr)
  glonbdr = glonbdr[wt2order]
  glatbdr = glatbdr[wt2order]
  lout2 = strcompress(string(round(glonbdr)), /rem)
  filename2 = 'grs-'+lout2+'-cube.fits'
  fn2 = grs_dir+filename2

  runspec = fltarr(1000)
  bdrspec = fltarr(1000)
  v_std = findgen(1000)*212/1e3-50.0


  for j = 0, n_elements(ind)-1 do begin 
     if not file_test(fn[j]) then continue
; Pick filename and load new if necessary
     if fn[j] ne lastfn then grs_data = readfits(fn[j], hd)
     lastfn = fn[j]
     extast, hd, astrom
     rdhd, hd, s = h
     sz_grs = size(grs_data)
     ad2xy, glon[j], glat[j], astrom, x, y
;; ; Look up spectrum
     if x lt 0 or x gt sz_grs[1]-1 or $
        y lt 0 or y gt sz_grs[2]-1 then begin
        wts[j] = 0
        continue
     endif
     spectrum = interpolate(grs_data,replicate(x,sz_grs[3]),replicate(y,sz_grs[3]),indgen(sz_grs[3]))

;; ; Interpolate to common velocity scale
;     v_std = findgen(1000)*h.cdelt[2]/1e3-50.0
     specout = interpol(spectrum, h.v, v_std)
     specout[where(v_std lt min(h.v) or v_std gt max(h.v))] = !values.f_nan
;     specout = smooth(specout,/edge_trun,11,/nan)
     runspec = runspec+specout*wts[j]
  endfor
  
  for j = 0,n_elements(ind2)-1 do begin
     if not file_test(fn2[j]) then continue
; Pick filename and load new if necessary
     if fn2[j] ne lastfn then grs_data = readfits(fn2[j], hd)
     lastfn = fn2[j]
     extast, hd, astrom
     rdhd, hd, s = h
     sz_grs = size(grs_data)
     ad2xy, glonbdr[j], glatbdr[j], astrom, x, y
;; ; Look up spectrum
     if x lt 0 or x gt sz_grs[1]-1 or $
        y lt 0 or y gt sz_grs[2]-1 then begin
        ctbdr = ctbdr - 1
        continue
     endif
     spectrum = interpolate(grs_data,replicate(x,sz_grs[3]),replicate(y,sz_grs[3]),indgen(sz_grs[3]))

     specout = interpol(spectrum, h.v, v_std)
     specout[where(v_std lt min(h.v) or v_std gt max(h.v))] = !values.f_nan
;     specout = smooth(specout,/edge_trun,11,/nan)

     bdrspec = bdrspec+specout
  endfor
  filter = savgol(32,32,0,6)

  specout = runspec/total(wts)
  spectrum = convol(specout,filter,/edge_trun,/nan) >0

  bdrspec = bdrspec/ctbdr
  bdrspec = convol(bdrspec,filter,/edge_trun,/nan) >0


; Final Spectrum     
  emmask = spectrum gt 4*mad(spectrum)
  spectrum = ((spectrum - bdrspec)>0)*emmask

;;   ps,/ps,file='specex2.ps',xsize=6.5,ysize=4.0,/def,/jour,$
;;      /color
;;   setcolors,/sys
;;   plot,v_std,spectrum,xrange=[-10,120],/xst,xtitle='V!DLSR!N (km/s)',$
;;        ytitle='!6T!DA!N (K)'
;; ;  oplot,v_std,bdrspec,ps=10,color=!red
;;   ps,/x  

  spectrum = spectrum^2


  if total(spectrum) eq 0 then return,prob

; Probability distribution for Kinetmatic distances.

  l = bgps.glon_max
  b = bgps.glat_max
  
  z = d*sin(b*!dtor)
  dproj = d*cos(b*!dtor)
  r = sqrt(r0^2+dproj^2-2*R0*dproj*cos(l*!dtor))

; Angle between GC line and R vector to object
  theta = acos((R0-dproj*cos(l*!dtor))/R)
; Angle between LOS and Velocity Vector
  phi = !pi/2-l*!dtor-theta

; This is VLSR for correct LSR
  vobs = (v0-vs)*cos(b*!dtor)*cos(phi)-v0*sin(l*!dtor)
  
  bigu = 10.0
  bigv = 5.23
  bigw = 7.17
  
  vhelio = vobs-(bigu*cos(l*!dtor)-bigv*sin(l*!dtor))*cos(b*!dtor)-$
           bigw*sin(b*!dtor)

  euler, l, b, ra, dec, 2
  solarmotion_ra = ((17+49/6d1+58.667/3.6d3)*15)
  solarmotion_dec = (28+7/6d1+3.96/3.6d3)
  solarmotion_mag = 16.55294
  gcirc, 2, solarmotion_ra, solarmotion_dec, ra, dec, theta2
; This is the VLSR that a telescope would report
  vlsr = vhelio+solarmotion_mag*cos(theta2/206265.)
  prob = dblarr(n_elements(vlsr)) 
  for q = 0,n_elements(v_std)-1 do begin 
     prob = prob + exp(-(vlsr-v_std[q])^2/(2*errvlsr^2))*spectrum[q]
  endfor
  prob = prob/total(prob)
;   if total(specout gt 0) gt 0 then stop
;    spec.grs_vel = v_std
;    spec.grs_spec = specout
  
  
  return,prob
end
